nt_detrend (original)

[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter,wsize) - robustly remove trend
 
  y: detrended data
  w: updated weights
  r: basis matrix used

  x: raw data
  order: order of polynomial or number of sin/cosine pairs
  w0: weights
  basis: 'polynomials' [default] or 'sinusoids', or user-provided matrix
  thresh: threshold for outliers [default: 3 sd]
  niter: number of iterations [default: 3]
  wsize: window size for local detrending [default: all]

 This NEW (circa Oct 2019) version of detrend allows detrending to be
 applied to smaller overlapping windows, which are then stitched together
 using overlap-add. This is not described in the paper.

0001 function [y,w,r]=nt_detrend(x,order,w0,basis,thresh,niter,wsize)
0002 %[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter,wsize) - robustly remove trend
0003 %
0004 %  y: detrended data
0005 %  w: updated weights
0006 %  r: basis matrix used
0007 %
0008 %  x: raw data
0009 %  order: order of polynomial or number of sin/cosine pairs
0010 %  w0: weights
0011 %  basis: 'polynomials' [default] or 'sinusoids', or user-provided matrix
0012 %  thresh: threshold for outliers [default: 3 sd]
0013 %  niter: number of iterations [default: 3]
0014 %  wsize: window size for local detrending [default: all]
0015 %
0016 % This NEW (circa Oct 2019) version of detrend allows detrending to be
0017 % applied to smaller overlapping windows, which are then stitched together
0018 % using overlap-add. This is not described in the paper.
0019 
0020 nt_greetings;
0021 
0022 %% arguments
0023 if nargin<2; error('!'); end
0024 if nargin<3; w0=[]; end
0025 if nargin<4||isempty(basis); basis='polynomials'; end
0026 if nargin<5||isempty(thresh); thresh=3; end
0027 if nargin<6||isempty(niter); niter=3; end
0028 if nargin<7; wsize=[]; end
0029 
0030 if iscell(x)
0031     if ~isempty(w0); error('not implemented'); end
0032     y={}; w={}; r={};
0033     for iTrial=1:numel(x);
0034         [y{iTrial},w{iTrial},r{iTrial}]=nt_detrend(x{iTrial},order,w0,basis,thresh,niter,wsize);
0035     end
0036     return
0037 end
0038 
0039 w=w0;
0040 if ~isempty(w)
0041     w=w(:);
0042     if numel(w)<size(x,1)
0043         % assume that w contains indices, set them to 1
0044         idx=w;
0045         w=zeros(size(x,1),1);
0046         w(idx)=1;
0047     end
0048 end
0049 
0050 if isempty(wsize) || ~wsize
0051     % standard detrending (trend fit to entire data)
0052     [y,w,r]=nt_detrend_helper(x,order,w,basis,thresh,niter);
0053 else
0054     wsize=2*floor(wsize/2);
0055     
0056     % do some sanity checks because many parameters
0057     if numel(order)>1; error('!'); end
0058     if ~isempty(w) && ~(size(w,1)==size(x,1)) ; disp(size(w)); error('!'); end
0059     if ~(isempty(basis) || isstring(basis) || ~(isnumeric(basis)&&size(basis,1)==size(x,1))); error('!'); end
0060     if thresh==0; error('thresh=0 is not what you want...'); end % common mistake
0061     if numel(thresh)>1; error('!'); end
0062     if numel(niter)>1; error('!'); end
0063 
0064     dims=size(x); nchans=dims(2);
0065     x=x(:,:); % concatenates dims >= 2
0066     w=w(:,:);
0067     if isempty(w); w=ones(size(x)); end
0068     if size(w,2)==1; w=repmat(w,1,size(x,2)); end
0069 
0070 
0071     % (1) divide into windows, (2) detrend each, (3) stitch together, (4)
0072     % estimate w
0073 
0074     for iIter=1:niter
0075 
0076         y=zeros(size(x));
0077         trend=zeros(size(x));
0078         a=zeros(size(x,1),1);
0079 
0080     %     figure(1); clf
0081 
0082         offset=0;
0083         while true
0084             start=offset+1;
0085             stop=min(size(x,1),offset+wsize);
0086 
0087             % if not enough valid samples grow window:
0088             counter=5;
0089             while any (sum(min(w(start:stop),2))) <wsize
0090                 if counter <= 0 ; break; end 
0091                 start=max(1,start-wsize/2);
0092                 stop=min(size(x,1),stop+wsize/2);
0093                 counter=counter-1;
0094             end
0095             if rem(stop-start+1,2)==1; stop=stop-1; end
0096             wsize2=stop-start+1;
0097 
0098             % detrend this window
0099             NITER=1;
0100             yy=nt_detrend_helper(x(start:stop,:),order,w(start:stop,:),basis,thresh,NITER);
0101 
0102             % triangular weighting
0103             if start==1
0104                 b=[ones(1,wsize2/2)*wsize2/2, wsize2/2:-1:1]';
0105             elseif stop==size(x,1)
0106                 b=[1:wsize2/2, ones(1,wsize2/2)*wsize2/2]';
0107             else
0108                 b=[1:wsize2/2, wsize2/2:-1:1]';
0109             end
0110 
0111             % overlap-add to output
0112             y(start:stop,:)=y(start:stop,:)+bsxfun(@times,yy,b);
0113             trend(start:stop,:)=trend(start:stop,:)+bsxfun(@times,x(start:stop,:)-yy,b);
0114 
0115             a(start:stop,1)=a(start:stop)+b;
0116 
0117             offset=offset+wsize/2;
0118             if offset>size(x,1)-wsize/2; break; end
0119         end
0120         
0121         if stop<size(x,1); y(end,:)=y(end-1,:); a(end,:)=a(end-1,:); end; % last sample can be missed
0122         
0123         y=bsxfun(@times,y,1./a); y(find(isnan(y)))=0;
0124         trend=bsxfun(@times,trend,1./a);  trend(find(isnan(trend)))=0;
0125 
0126         % find outliers
0127         d=x-trend; 
0128 
0129 
0130         if ~isempty(w); d=bsxfun(@times,d,w); end
0131         ww=ones(size(x));
0132         ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0133         clear d
0134 
0135         % update weights
0136         if isempty(w); 
0137             w=ww;
0138         else
0139             w=min(w,ww);
0140         end
0141         clear ww;
0142 
0143     end % for iIter...
0144     
0145     w=[];r=[]; % not informative
0146     
0147 end % if isempty(wsize)...
0148 
0149 if ~nargout
0150     % don't return, just plot
0151     subplot 411; plot(x); title('raw'); xlim([1,size(x,1)])
0152     subplot 412; plot(y); title('detrended'); xlim([1,size(x,1)])
0153     subplot 413; plot(x-y); title('trend'); xlim([1,size(x,1)])
0154     subplot 414; nt_imagescc(w'); title('weight');
0155     clear y w r
0156 end
0157 end
0158 
0159 %% Original version of detrend (no windows) is called by new version (windows)
0160 function [y,w,r]=nt_detrend_helper(x,order,w,basis,thresh,niter)
0161 %[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter) - robustly remove trend
0162 %
0163 %  y: detrended data
0164 %  w: updated weights
0165 %  r: basis matrix used
0166 %
0167 %  x: raw data
0168 %  order: order of polynomial or number of sin/cosine pairs
0169 %  w: weights
0170 %  basis: 'polynomials' [default] or 'sinusoids', or user-provided matrix
0171 %  thresh: threshold for outliers [default: 3 sd]
0172 %  niter: number of iterations [default: 3]
0173 %
0174 % Noise tools
0175 % See nt_regw().
0176 %
0177 % The data are fit to the basis using weighted least squares. The weight is
0178 % updated by setting samples for which the residual is greater than 'thresh'
0179 % times its std to zero, and the fit is repeated at most 'niter'-1 times.
0180 %
0181 % The choice of order (and basis) determines what complexity of the trend
0182 % that can be removed.  It may be useful to first detrend with a low order
0183 % to avoid fitting outliers, and then increase the order.
0184 %
0185 % Examples:
0186 % Fit linear trend, ignoring samples > 3*sd from it, and remove:
0187 %   y=nt_detrend(x,1);
0188 % Fit/remove polynomial order=5 with initial weighting w, threshold = 4*sd:
0189 %   y=nt_detrend(x,5,w,[],4);
0190 % Fit/remove linear then 3rd order polynomial:
0191 %   [y,w]=nt_detrend(x,1);
0192 %   [yy,ww]=nt_detrend(y,3);
0193 %
0194 nt_greetings;
0195 
0196 % arguments
0197 if nargin<2; error('!'); end
0198 if nargin<3; w=[]; end
0199 if nargin<4||isempty(basis); basis='polynomials'; end
0200 if nargin<5||isempty(thresh); thresh=3; end
0201 if nargin<6||isempty(niter); niter=3; end
0202 
0203 if thresh==0; error('thresh=0 is not what you want...'); end % common mistake
0204 
0205 dims=size(x);
0206 x=x(:,:); % concatenates dims >= 2
0207 w=w(:,:);
0208 
0209 if size(w,2)==1; w=repmat(w,1,size(x,2)); end
0210 
0211 %% regressors
0212 if isnumeric(basis)
0213     r=basis;
0214 else
0215     switch basis
0216         case 'polynomials'
0217             r=zeros(size(x,1),numel(order));
0218             lin=linspace(-1,1,size(x,1));
0219             for k=1:order
0220                 r(:,k)=lin.^k;
0221             end
0222         case 'sinusoids'
0223             r=zeros(size(x,1),numel(order)*2);
0224             lin=linspace(-1,1,size(x,1));
0225             for k=1:order
0226                 r(:,2*k-1)=sin(2*pi*k*lin/2);
0227                 r(:,2*k)=cos(2*pi*k*lin/2);
0228             end
0229         otherwise
0230             error('!');
0231     end
0232 end
0233 
0234 
0235 % remove trends
0236 % The tricky bit is to ensure that weighted means are removed before
0237 % calculating the regression (see nt_regw).
0238 
0239 for iIter=1:niter
0240     
0241     % weighted regression on basis
0242     [~,y]=nt_regw(x,r,w);
0243     
0244     % find outliers
0245     d=x-y; 
0246     if ~isempty(w); d=bsxfun(@times,d,w); end
0247     ww=ones(size(x));
0248     ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0249      
0250     % update weights
0251     if isempty(w); 
0252         w=ww;
0253     else
0254         w=min(w,ww);
0255     end
0256     clear ww;    
0257 end
0258 y=x-y;
0259 y=reshape(y,dims);
0260 w=reshape(w,dims);
0261 
0262 end
0263 
0264 %% test code
0265 function test_me
0266 if 0
0267     % basic
0268     x=(1:100)'; x=x+ randn(size(x));
0269     WSIZE=30;
0270     y=nt_detrend2(x,1,[],[],[],[],WSIZE);
0271     figure(1); clf; plot([x,y]);
0272 end
0273 if 0
0274     % basic
0275     x=(1:100)'; x=x+ randn(size(x));
0276     x(40:50)=0;
0277     WSIZE=30;
0278     [yy,ww]=nt_detrend2(x,1,[],[],2,[],WSIZE);
0279     [y,w]=nt_detrend(x,1);
0280     figure(1); clf; subplot 211; 
0281     plot([x,y,yy]);
0282     subplot 212; plot([w,ww],'.-');
0283 end
0284 if 0
0285     % detrend biased random walk
0286     x=cumsum(randn(1000,1)+0.1);
0287     WSIZE=200;
0288     [y1,w1]=nt_detrend(x,1,[]);
0289     [y2,w2]=nt_detrend2(x,1,[],[],[],[],WSIZE);
0290     figure(1); clf; 
0291     plot([x,y1,y2]); legend('before', 'after');
0292 end
0293 if 0
0294     % weights
0295     x=linspace(0,100,1000)';
0296     x=x+3*randn(size(x));
0297     x(1:100,:)=100;
0298     w=ones(size(x)); w(1:100,:)=0;
0299     y=nt_detrend2(x,3,[],[],[],[],WSIZE);
0300     yy=nt_detrend2(x,3,w,[],[],[],WSIZE);
0301     figure(1); clf; plot([x,y,yy]); legend('before', 'unweighted','weighted');
0302 end
0303 if 0
0304     [p,x]=nt_read_data('/data/meg/theoldmanandthesea/eeg/mc/MC_aespa_speech_43.mat'); x=x'; x=x(:,1:128); %x=x(1:10000,:);
0305     %[p,x]=nt_read_data('/data/meg/arzounian/ADC_DA_140521_p20/ADC_DA_140521_p20_01_calib'); x=x'; x=x(1:10000,:);
0306     
0307     x=nt_demean(x);
0308     figure(1);
0309     nt_detrend(x,3);
0310     figure(2);
0311     WSIZE=1000;
0312     nt_detrend2(x(:,:),3,[],[],[],[],WSIZE);
0313 end
0314 end
0315 
0316